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Abstract. Markov chains are a convenient means of generating real- 
izations of networks, since they require little more than a procedure for 
rewiring edges. If a rewiring procedure exists for generating new graphs 
with specified statistical properties, then a Markov chain sampler can 
generate an ensemble of graphs with prescribed characteristics. However, 
successive graphs in a Markov chain cannot be used when one desires in- 
dependent draws from the distribution of graphs; the realizations are 
correlated. Consequently, one runs a Markov chain for N iterations be- 
fore accepting the realization as an independent sample. In this work, we 
devise two methods for calculating N. They are both based on the binary 
"time-series" denoting the occurrence/non-occurrence of edge (u,v) be- 
tween vertices u and v in the Markov chain of graphs generated by the 
sampler. They differ in their underlying assumptions. We test them on 
the generation of graphs with a prescribed joint degree distribution. We 
find the N oc \E\, where \E\ is the number of edges in the graph. The two 
methods are compared by sampling on real, sparse graphs with 10 3 — 10 4 
vertices. 

Keywords: graph generation, Markov chain Monte Carlo, independent 
samples 

1 Introduction 

Markov chain Monte Carlo (MCMC) methods are a common means of generat- 
ing realizations of graphs which share similar characteristics since they require 
nothing more than a procedure that can generate a new graph by "rewiring" 
the edges of an existing graph. Much of their use to date has been in gener- 
ating graphs with a prescribed degree distribution 1 2 3 4J. Other efforts have 
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used MCMC to generate graphs with a prescribed joint degree distribution [5]. 
MCMC methods require a graph to start the chain; thereafter, the "rewiring" 
procedure generates new, realizations which preserve certain graph character- 
istics. The specific characteristic (s) that are preserved depend entirely on the 
"rewiring" procedure. 

MCMC methods for generating graphs has two drawbacks. Initialization bias 
arises from the fact that the starting graph may not even lie in the population 
of graphs that we seek to sample, or may be an outlier in that population. The 
second issue, autocorrelation in equilibrium, arises from the fact that successive 
samples drawn by the MCMC sampler are correlated and an empirical distribu- 
tion constructed from them would result in the statistical error (variance) to be 
2ri„tlarger than the distribution constructed using independent samples. Here 
T int , the integrated autocorrelation time, is a measure of how slowly correlation 
in graphical metrics, calculated from an MCMC series, decays; see Sections 2 
and 3 in Sokal's lecture notes [BJ. 

Sokal's method [BJ for deciding the "sufficiency" of samples obtained from 
MCMC revolve around autocorrelation. The method is general, and was adapted 
for use with graphs in [5]. Consider an edge (it, v) between labeled vertices u 
and v in the ensemble of graphs generated by the MCMC chain. Denoting its 
occurrence/non-occurrence in the chain of graphs by 1/0 gives us a binary time- 
series {Z t }, t = 1 . . . T, with an empirical mean fj,. The auto-correlation, with lag 
I, is given by C(l) — ({Z t } — p)({Z t +i} — fx),t = 1...T — I and the normalized 
version of it, p(l) = C(l)/C(0) can be used to gauge whether the autocorrelation 
in the time-series is observed to be decreasing. In [5] , the authors used this metric, 
applied to all edges in the graphs that were sampled, to ensure that their MCMC 
chain was mixed. One can also set, loosely speaking (for details, consult [BJ), a 
minimum threshold p m im identify the corresponding lag l m in, and retain every 
l^in entry in the MCMC chain to serve as independent samples. However, this 
method has two practical drawbacks. First the autocorrelation analysis has to 
be performed for all the edges (potentially, \V\ 2 in number) that might appear 
in the MCMC chain, which quickly becomes prohibitively expensive for large 
graphs. Secondly, it requires a user input, p m im which may have an arbitrary 
effect on graphical properties of the ensemble. These shortcomings motivate our 
work. 

In this paper, we propose two different methods for generating independent 
graphs using an MCMC method. The first, which we call Method A or "multiple 
short runs" , determines the number of iterations N an MCMC method has to be 
run to "forget" the initial graph and minimize the initialization bias. The second 
approach, Method B or "one long run" , requires K MCMC iterations. This long 
run is thinned by a factor k (i.e., every k th MCMC iteration is preserved) to gen- 
erate K/k independent samples. Both methods are intended to be approximate, 
but simple to evaluate, so that they can be employed in practice to gauge the 
"sufficiency" of MCMC iterations. The two methods for extracting independent 
graphs are tested on an MCMC chain with the setup described in [5] . We explore 
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the practical impact of approximations in our methods. We restrict ourselves to 
undirected graphs with labeled nodes. 

In the next section (Sec. [2]) we describe the procedure used to "rewire" a 
graph to create a new graph realization with the same joint degree distribution. 
In Sec. [3] we describe the two methods for generating independent samples. In 
Sec. |1J we test the methods on real sparse graphs. We conclude in Sec. [5j 

2 A Markov chain algorithm for sampling graphs with a 
given joint degree distribution 

Consider an undirected graph G — (V,E), where |V| = n and \E\ — m. The 
degree distribution of the graph is given by the vector /, where f(d) is the 
number of vertices of degree d. The joint degree distribution list the number of 
edges incident between vertices of specified degrees. Formally, the n x n matrix 
J denotes the joint degree distribution, where the entry J(i,j) is the number of 
edges between vertices of degree i and degree j. Stanton and Pinar [5] studied 
the problem of generating and random sampling a graph with a given joint 
degree distribution. They proposed a greedy algorithm to construct an instance 
of a graph with a specified (feasible) degree distribution, as well as a Markov 
chain algorithm to generate random samples of graphs with the same degree 
distribution. 

For the purposes of this paper, we will only focus on the Markov chain algo- 
rithm. The rewiring operation that moves us between the nodes of the Markov 
chain is depicted in Fig. [lj At the first step, one picks an edge (ui,v) at ran- 
dom and thereafter, one of the end vertices, e.g., u\. We wish to break (u\,v) 
and connect u\ and v to others without violating the prescribed joint degree 
distribution. We, therefore, search for another edge (u2,w) where d U2 = d Ul or 
d w = d Ull where d p denotes the degree of node p. WLOG, let d U2 — d Ul . Swap- 
ping the edges i.e. creating edges (ui,w) and (u2,v) while destroying (ui,v) 
and (u2,w) leaves the joint degree distribution unchanged while changing the 
connectivity pattern of the graph. If the resulting graph is simple, the graph is 
retained by the MCMC chain. 

The procedure results in inter-edge correlation. A particular edge is chosen for 
swapping, on average, once every \E\ Markov chain iterations. In [5], this proce- 
dure was used to generate graphs (of moderate sizes, with \V\ < 23, 000), using an 
MCMC sampler. Autocorrelation analysis showed that the Markov chain mixes 
and the autocorrelation decays with edge-dependent rates. Empirically it was 
observed that an edge (u,v) de-correlated slowly if J(d u ,d v )/(fd u fd v ) ~ 0.5. 
Empirically, it was observed that the autocorrelations of the edges decreased 
very sharply. 

3 Methods for calculating independence 

The algorithm in Sec. [2] has two variants - one where the resulting graph may 
not be simple, and another where the graph was always simple (A graph is 



U 1 u 1 u 2 



Step 1: Pick an edge 
(u 1; v), and pick one 
of its vertices, e.g., u 1 



Step 2: Pick another edge 
(u,w), such that d(u 1 )=d(u 2 j 
or d(uj=d(w) 




Step 3: Swap edges 



Fig. 1. The swapping operation for the Markov chain algorithm. 



simple if it does not have self-loops or parallel edges). For the first variant, 
earlier work [T] implies that this chain has a polynomial mixing time. The second 
variant is not even known to have any bounds on the mixing time. As discussed 
in Sec.[TJ one would like to provide a length to run the MCMC that, although not 
guaranteeing complete mixing, at least gives some confidence that the sampled 
graph is fairly "random." To that end, we will we will approximate the behavior 
of a single edge by a Markov chain. We stress that we do not give a proof, but 
only a mathematical argument justifying this. Below we present two methods 
for generating independent graph samples. 



3.1 Method A - "multiple short runs" 

We refer to this method as Method A or the "multiple short runs" method. 
We generate M graph samples by running M independent Markov chain for N 
iterations before accepting the resulting graph. All the chains are run from the 
same initial graph; however, the state in the random number generator in each 
of the M MCMC chains are distinct. 

Consider a fixed pair of labeled vertices {u,v}. We will approximate the 
occurrence of an edge between (u, v) as a two state Markov chain. Note that this 
is an approximation, since these transitions depend on the remaining graph if 
J is to be preserved. Nonetheless, these dependencies appear to be weak. The 
first coordinate of the matrix is state (no edge) and the second coordinate is 
1, indicating the existence of an edge. The transition matrix T for this chain is 

rp _ I 1 — a iJ a i,j \ I'll 

{ fa i W' [ ) 

where i = d u and j = d v are the degrees of vertices u and v. a i : j and fiij are 
positive fractions and (ctij + ftj) < 1. The eigenvalues of the transition matrix 
are 1 and 1 — (ctjj + Below, we construct a model for ctjj and /3ij. 

Suppose the state is currently 0. The state will become 1 if the edge (u, v) is 
swapped in. Let the two edges chosen by the algorithm be e and e', in that order. 
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The edge (u,v) is swapped in if e contains u and e' contains v (or vice versa). 
Furthermore, the endpoint u must be chosen, and the other end of e' must have 
degree d u . The probability that e contains u and u is chosen as an endpoint 
is exactly d u /2m. The probability we choose the edge e' that is incident to v 
depends on the number of neighbors of v whose degree is d u . Clearly, this depends 
on the graph structure (leading in a non-Markov probability of this transition). 
We heuristically guess this number based on the joint degree distribution. The 
number of edges from degree d u to degree d v vertices is J(d u , d v ). Of these, the 
average number of edges incident to a fixed vertex of degree d v is 3(d Ul d v )/f v . 
We shall approximate the number of edges incident to v with the other endpoint 
of degree d u by this quantity. The total probability is 

du ^ 3{du:d v ) d u J(d u ,d v ^) 

2m mf(d v ) 2m 2 f(d v ) 

The edge (u, v) is also swapped in when the reverse happens (so we choose v as 
an endpoint, and an edge incident to u with the other endpoint of degree d v ). 
The total transition probability from to 1 is approximated by 



d u J(d U7 d v ^ d v J(d u ,d v ^ J(d ii ,d l; ) / d u d v 
* iJ = 2m 2 f(d v ) 2m 2 f(d u ) = 2m 2 [j(d^ J{du) 



(2) 



We now address the transition from 1 to 0. Suppose (u, v) is currently an edge. 
If the first edge e is chosen to be (u, v), then (u,v) will definitely be swapped 
out. The probability of this is 1/m. If the random endpoint chosen has degree 
d u (and is not u), then we might choose e' to be (u,v). The total probability of 
this is 

(/K) - IK x 1 = (f(du) - l)d u 
2m m 2m? 

The roles of u and v can also be reversed, so the total transition probability from 
1 to is 

f(d. u )d u + f(d v )d v - d u - d v 
2m 2 

and so 

J_ f(d u )d u + f{d v ) d v 
m 2m 2 
We proceed to determining the number of iterations N to run the Markov 
chain. We start the Markov chain Ai with an initial distribution v (which 
is either (0,1) or (1,0)). Ai, which is represented by Ti.j (Eq. [T|, is run for 
N = ln{\ I e) I (a + fi) iterations, e > 0. We have dropped the subscripts i and j, 
since it is implied that this model is being derived for an edge (it, v) with vertices 
of degrees i and j. After N steps, we realize a 2-state distribution p = (po,Pi), 
which is different from the stationary distribution be u = (uq,Ui). 

Denote the unit 2-norm eigenvectors of T, corresponding to the eigenvalues 1 
and 1 — (a+/3), as ei and e2. The initial state can be expressed as v = ciei+C2e2. 
After N applications of the transition matrix we get 

p = T w v = Cl T N ei + c 2 T N e 2 = ciei + c 2 (1 - (a + /3)) W e 2 . 



() 



Since (1 — {a+/3}) < 1, the second term decays with N and ciei is the stationary 
distribution. We can bound the decaying term as 

||(1 - (a + /3)) JV c 2 e 2 || 2 = (1 - (a + /3)) ln(1/£)/(a+ «c 2 ||e 2 || 2 < exp(- ln(l/e)) = e 

Hence, ||p — u|| 2 < e, and so each \pi — Ui\ is at most e. Further, from Eq.[2]and 
[3j we see that a + > 1 /m (to leading order) and consequently 

N = ln ^f < mln(l/e) = \E\ In (I) . (4) 
3.2 Method B - "one long run" 

We propose a second method, which we refer to as Method B or "one long run", 
for generating independent graphs. The procedure involves running a Markov 
chain for a large number of steps K and thinning it by a factor k i.e., pr eserv ing 
every k th instance of the chain. Comparing with the development in Sec. 3.1 we 
expect k ~ N . 



Similar to Method A (Sec. |3.1[ ), this method too begins with the binary 
time-series of edge occurrence {Z t }. As observed in [5], the autocorrelation in 
{Z t } decays for all edges. Consequently it is possible to successively thin the 
chain {Z t } (i.e., retain every k th element to obtain {Z^}, the k— thinned chain) 
and compare the likelihoods that the chains were generated by (1) independent 
sampling or (2) by a first-order Markov process. When sufficiently thinned, the 
independent sampling model is expected to fit the data better. Using this as 
the stopping criterion removes an ambiguity (user-specified tolerances). We will 
employ a method based on comparison of log-likelihoods of model fit. We de- 
rive these expressions below. While this technique has been applied in other 
domains |7l8j . but this paper is the first application of this technique to graphs. 

Consider the chain {Z^}. We count the number, Xij, of the £ (0,1) 

transitions in it. are used to populate X, a 2 x 2 contingency table. Dividing 
each entry by the length of thinned chain K/ k—1 provides us with the empirical 
probabilities pij of observing an transition in {Z^}. Let pTj and x~ij = 

(K/k — l)pTj be the predictions of the probabilities and expected values of the 
table entries provided by a model. In such a case, the goodness-of-fit of the model 
is provided by a likelihood ratio statistic (called the G 2 -statistic; Chapter 4.2 
in [5]) and a Bayesian Information Criterion (BIC) score 

° 2 = - 2 E E *H l0 S (— V BIC = ° 2 + n l0 § ( T J > ) 

where n is the number of parameters in the model used to fit the table data. 
Typically log-linear models are used for the purpose (Chapter 2.2.3 in [S]); the 
log-linear models for table entries generated by independent sampling and a 
first-order Markov process are 

(6) 
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where superscripts /, M indicate an independent and Markov process respec- 
tively. The maximum likelihood estimates (MLE) of the model parameters (w^x) 
are available in closed form (Chapter 3.1.1 in They lead to the model pre- 
dictions below 



and xf' = Xij, (7) 



a x++ — -*v> 

where Xi + and x + j are the sums of the table entries in row i and column j 
respectively. x++ is the sum of all entries (i.e., K/k—1, the number of transitions 
observed in {Z^}, or the total number of data points). We compare the fits of 
the two models thus: 



i=l i=l 



ABIC = BICW - BICM = x » log I ?T I ~ log ( T ~ 1 ) ' (8) 

= 4=0 




c\j = Xij and the fact that the log-linear model 
for a Markov process has one more parameter than the independent sampler 
model. Large BIC values indicate a bad fit. A negative ABIC indicates that an 
independent model fits better than a Markov model. 

The procedure for identifying a suitable thinning factor k then reduces to 
progressively thinning {Z^} till ABIC in Eq. [8] becomes negative. We search 
for k in powers of 2. The value of k so obtained varies between edges and con- 
servatively, we take the largest k, fc*. However, this may be too conservative, i.e, 
fc* 3> N, if a few edges are seen to display a slow autocorrelation decay. If we 
are interested in certain global metrics for graph e.g., maximum eigenvalue etc, 
a few correlated edges are unlikely to have any substantial effect. Thus, one may 
be able to thin with a k ~ N <C fc*. We will test this empirically in Sec. [4] 



4 Tests with real graphs 



In this section we first explore the impact of e (as defined in Sec. 3.1) on the 



ensemble of graphs generated by a Markov chain, and choose a e for further use. 



Thereafter we compare the graphs generated by the Methods A and B (Sec. 3.1 
and Sec. |3.2[ ) and gauge the impact of choosing a thinning factor k < k*. All 
tests are done with four real networks - the neural network of C. Elegans [TO] 
(referred to as "C. Elegans"), the power grid of the Western states of US [TU] 
(called "Power"), co-authorship graph of network science researchers (TTj (re- 
ferred to as "Netscience" ) and a 75,000 vertex graph of the social network at 
Epinions.com [12] ( "Epinions" ) . Their details are in Table[TJ The first three were 
obtained from [13] while the fourth was downloaded from [M]. All the graphs 
were converted to undirected graphs by symmetrizing the edges. 

We start the Markov chain using the real networks listed in Table [T] When 
comparing ensembles of graphs, we will use the (distributions of) global cluster- 
ing coefficient, number of triangles in the graphs, the graph diameter and the 
maximum eigenvalue as metrics. 



Table 1. Characteristics of the graphs used in this paper. (| V\, \E\) are the numbers 
of vertices and edges in the graph, N are the number of Markov chain steps used for 
generating graphs in Sec. |3.1| and k is its equivalent obtained by the method in Sec. |3.2| 
K/k* are the number of graph samples, obtained by thinning a long run, that were 
used to generate distributions in the figures. 



Graph name 


(\V\,\E\) 


N/\E\ 


k,/\E\ 


K/k, 


C. Elegans 


(297, 4296) 


10 


13 


3582 


Netscience 


(1461, 5484) 


10 


49 


737 


Power 


(4941, 13188) 


10 


13 


1214 


Epinions 


(75879, 405740) 


30 


720 


various 



In Fig. [2] we investigate the impact of e in Method A ("many short runs"). 
We generate 1000 samples by running the Markov chain for 1\E\, 5|2?|, 1Q\E\ and 
15|-E| Markov chain iterations, corresponding to e = 0.37,6.7 x 10~ 3 ,4.5 x 10 -5 
and 3.06 x 10~ 7 . In Fig. [2J we plot the distributions for the first three graphs (in 
Table [l| and find that for all three, e < 5.0 x 10~ 3 lead to distributions which 
are very close. We will proceed with e = 4.5 x 10 -5 i.e., when we use Method A, 
we will mix the Markov chain 10|i?| times before extracting a sample. 

In Table [l] we see that Method B ( "one long run" ) method often prescribes a 
thinning factor that is larger than the one obtained using Method A ( "multiple 
short runs"). This large number is often due to the lack of autocorrelation decay 
in a few edges. We investigate whether such a lack has a significant impact on 
the graphical metrics that we have chosen. In Fig. [3] we plot distributions of the 
same metrics for the three graphs. The thinning factors are in Table [T] We see 
that the distributions are close, i.e., the existence of a few edges whose time- 
series are still correlated do not impact the metrics of choice. We have repeated 
these tests with other metrics and the same result holds true. 

We now address a large graph (Epinions). Since potentially \V\ 2 distinct 
edges might be realized during a Markov chain, it is infeasible to calculate a 
thinning factor for all the edges. Consequently, we perform the thinning analysis 
for only 0.1|i?| (40,574) edges, chosen randomly from all the distinct edges that 
are realized by the Markov chain. In Fig. [3] we plot the distribution of k obtained 
from the 40,574 sampled edges. We see that most of the k lie between 10|-E| and 
100|£|; edges with thinning factors outside that range are about two orders of 
magnitude less abundant. It is quite conceivable that there are edges (which 
were not captured by the sample) that would prescribe an even higher thinning 
factor. In order to check whether these edges have a significant impact on the 
distribution of graph metrics, we check their convergence as a function of the 
thinning factor. 

We generate separate ensembles of graphs. The reference ensemble is gener- 
ated using Method A, with N = 30|_E|. As seen in Table [l] certain edges will 
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Fig. 2. Plots of the distributions of the global clustering coefficient, the number of 
triangles in the graphs, the graph diameter and the max eigenvalue of the graph Lapla- 
cian for "Netscience" (left) and "Power" (right), evaluated after 1\E\, 5\E\, W\E\ and 
15\E\ iterations of the Markov chain (green, blue, black and red lines respectively). 
The corresponding values of e are in the legend. We see that the distributions converge 
at e < 1.0 _s . 



still be correlated (N = 720\E\ would make them independent). We then use 
Method B to generate graph ensembles with thinning factors k/\E\ < 720 which 
are also multiples of N. In Fig. [4] (right), the diameter distribution obtained with 
Method A is compared to that obtained with Method B. While the distributions 
are very similar, they do display some small differences. This is surprising since 
N = 30\E\ indicates a minuscule e. In addition, distributions obtained with 
k = 5N, 9N, and 137V show some differences between themselves, indicating 
that the edges that have not become independent have a small, but measurable 
impact on the graph diameter. Further, the distributions using smaller values of 
k are marginally wider (have a larger variance), indicating that they were con- 
structed using samples which were not completely independent. However, the 
differences are minute, and for practical purposes the graph ensemble generated 
using Method A with N = 30\E\ is identical to the one generated using Method 
B, per our chosen metrics. Consequently, despite its approximations, the results 
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Fig. 3. Comparison of the distributions of the global clustering coefficient, the number 
of triangles in the graphs, the graph diameter and the max eigenvalue of the graph 
Laplacian for "C. Elegans" (left), "Netscience" (middle) and "Power" (right), evalu- 
ated using Methods A ("many short runs") and B ("one long run"). We see that the 
distributions are very similar. The kernel density estimation used to generate the dis- 
tributions sometimes causes nonsensical artifacts e.g., a small, but negative clustering 
coefficient. For Method A, the Markov chain was run for 10\E\ iterations. Thinning 
factors for Method B are in Tabled] 



in Sec. 3.1 furnish a workable estimate of N, if one uses e < 10~ 5 . Further, fc* 
is generally too conservative if our aim is to obtain "converged" distributions of 
certain graph metrics. This arises from a few edges that de-correlate slowly, but 
have little effect on global graphical metrics due to their rarity. 



5 Conclusions 

We have developed a method that allows one to generate a set of independent 
realizations of graphs with a prescribed joint degree distribution. The graphs are 
generated using an MCMC approach, employing the algorithm described in [5] 
as the "rewiring" mechanism. The graphs so generated are tightly correlated; 
our two methods address the question of how one can decorrelate the chain. 
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Fig. 4. Left: The normalized thinning factor k/\E\ for the Epinions graph, as calcu- 
lated for the 40,574 sampled edges. We see that the most thinning factors are lie in 
(10|-E|, 100|.E|). Right: Plot of the graph diameter and distribution generated using 
Method A (with N = 30\E\) and Method B (with k equal to various multiples of N). 
We see that the distributions are very similar. 



The first method, variously called Method A or "multiple short chains", 
involves running the Markov chain for N steps before extracting a graph real- 
ization; the Markov chain is run repeatedly to generate samples. We developed 
a model (and a closed-form expression) to estimate N that allows the Markov 
chain to converge to its stationary distribution before a graph realization is ex- 
tracted from it. This model assumes that edges are independent. In reality, their 
behavior is correlated, which leads us to incur small errors. 

The second method, variously called Method B or "one long chain", is a 
data driven method. It uses the time-series of the occurence/non-occurence of 
edges in an MCMC run. It does not assume a constant joint degree distribution. 
It progressively thins the time-series (by retaining every k th element) and fits a 
first-order Markov and an independent sampling model to the data. The thinning 
process stops when the independent model has a higher likelihood (strictly, a 
lower BIC score) than the Markov process. Since this method is data-driven and 
does not require any user-defined tolerances, we use it to validate Method A. 
The method is not new, but does not seem to have been used in the generation 
of independent graphs. 

Comparing the two methods, we find that for practical purposes, the ensem- 
bles generated using Method A are statistically similar to those obtained with 
Method B, as gauged by a set of graph metrics. Even at tight tolerance values, 
a small number of edges in the graphs generated by Method A remain corre- 
lated, and the metrics' distributions are slightly wider. This problem is very 
small (nearly unmeasureable) in small graphs, but becomes measureable, but 
still small, for large graphs. 

While this work enables the generation of independent graphs, including 
large ones, it poses a number of questions for further investigation. For example, 



12 



being able to estimate or bound the difference in the distributions generated by 
Methods A and B would be helpful. Further, an intelligent way of identifying 
hard-to-decorrelate edges would reduce the computational burden of checking 
for the stopping criterion using Method B; currently, we simply use a random 
set of edges. Finally, it would be interesting if Method A could be extended to 
the generation of independent graphs when some graph property, other than the 
joint degree distribution, is held constant. This is currently being studied. 
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